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A FOUR-CYLINDER STIRLING ENGINE COMPUTER PROGRAM 
WITH DYNAMIC ENERGY EQUATIONS 


Carl J. Daniele and Carl F. Lorenzo 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


SUMMARY 


A computer program for simulating a four-cylinder Stirling engine is 
presented. The thermodynamic model includes both continuity and energy 
equations and simplified, first-order momentum terms (flow resistance). 

Drive dynamics and vehicle load effects are included. The computer program 
includes a model of a hydrogen supply system to accelerate the engine. 

The simulation can generate, at the user's option, steady-state or tran- 
sient data. Steady-state power and torque predictions compare well with ex- 
perimental data over a wide range of engine speeds and pressures. Simulated 
acceleration data compare well with results from a simplified four-cylinder 
model that was previously developed by the authors. Complete flow charts 
are provided. 


INTRODUCTION 


This report documents a computer program for simulating the steady-state 
and transient performance of a Stirling engine and serves as a users manual 
for the computer program. The evolution of the model is described below. 

Stirling engine simulations have been developed for predicting steady- 
state engine performance. The simulations are based on the solution of the 
governing aerothermodynamic equations and are therefore quite complex. 
Pressure drops due to fluid flow resistance within the engine and to heat 
transfer coefficients between the heat exchangers and fluid and between 
fluid and the casing are calculated continuously during a cycle. Stirling 
engine performance is very sensitive to these factors. 

Usual Stirling engine models include only two pistons and the working 
space between them. For analysis, the working space is usually segmented 
into control volumes corresponding to engine components - expansion space, 
heater, regenerator, cooler, and compression space. More than one control 
volume may be used to describe each component. Also, various heat transfer 
paths are included in the model. While many Stirling engines have more than 
one working space, the models used to predict performance generally do so 
for a single working space, and the resultant output power is multiplied by 
the number of working spaces. One such model has been developed by Tew, 
Jeffries, and Miao (ref. 1). It contains 13 control volumes in the working 
space. In that model, the energy equation is integrated within each control 
volume, while the pressure is assumed constant throughout the working space 



over a cycle. Pressure drops throughout the engine are then back-calculated 
once the temperature distribution is known. 

Although this approach is acceptable for performance calculations, it 
has shortcomings when the constant pressure assumption cannot be used, such 
as when the amount of hydrogen within a working space varies during a cycle. 
In a multicylinder engine, this can occur when transients are run from one 
operating point to another. Daniele and Lorenzo (ref. 2), as a first step 
in deriving a transient engine model, developed a model that can account for 
a change in the hydrogen-stored mass during a cycle. As in the Tew model, 
only one working space and two pistons are simulated. However, the number 
of control volumes in the working space is reduced to seven, one each for 
the expansion space, the heater, the cooler, and the compression space, and 
three for the regenerator. In that model, a thermal time constant associ- 
ated with the regenerator mesh is included. Flow resistances and heat 
transfer coefficients are held constant over an engine cycle. Values of 
these parameters are specified as input data, using average values over a 
cycle, calculated from the output of lew's model. Within each control vol- 
ume, the continuity and energy equations are integrated and a simplified, 
first-order momentum term (flow resistance) is calculated. Upwind differen- 
cing (ref. 3) is used to calculate the interface volume temperatures for use 
in the energy equation. The resultant model has 17 state variables and uses 
a backward-difference integration scheme for problem solution. Results gen- 
erated by that model compare well with experimental power and torque data 
over a wide range of speeds and pressures. 

The aforementioned modeling approach (i.e., those in which only one 
working space is simulated) produces acceptable steady-state (power and 
torque) predictions when the working spaces in the engine are isolated. 
However, proposed controls schemes for four-cylinder Stirling engines re- 
quire fluid transfer between working spaces. A more complete four-cylinder 
model is needed to analyze these control schemes. 

One such four-cylinder model is also described in reference 2. In that 
model, each working space is further simplified to reduce computer run time. 
This simplification involves (1) reducing the number of control volumes in a 
working space from seven to three, (2) keeping the temperature within a con- 
trol volume constant, and (3) assuming that the gas and regenerator mesh 
temperatures are equal. Steady-state results from this model are presented 
in references 2 and 4 and agree reasonably well with experimental data. The 
model was used to study various transient phenomena (ref. 4), and the asso- 
ciated computer program is documented in reference 5. 

The simplified four-cylinder model is well suited for controls studies 
since it can produce transient results reasonably fast. Thus, many differ- 
ent cases can be run without using excessive computer time. However, the 
model lacks temperature variability, which may be important when simulating 
proposed control schemes, such as the injection of working fluid at one tem- 
perature into a control volume with resident working fluid at a different 
temperature. 

To overcome this deficiency, the seven-volume, single-working-space 
model discussed earlier was extended to include all four pistons and working 
spaces; the resultant model is the subject of this report. A backward- 
difference integration scheme was used to account for the widely varying 
time constants associated with the flow and temperature dynamics. The inte- 
gration scheme makes use of a multivariable Newton-Raphson iteration method 
(ref. 6) for convergence at a time point. To speed calculations, the Broy- 
den update algorithm (ref. 7) was used to update the inverted Jacobian ma- 
trix during a transient. 
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This paper documents the seven-volume, four-working-space model and its 
computer implementation. A users manual and a test case are provided, along 
with complete flow charts. Results from the simulation are compared with 
results from the simplified four-cylinder model. 


MODEL DESCRIPTION 


The model consists of three parts: thermodynamic, drive geometry, and 

hydrogen supply. A schematic of the seven-volume, four-cylinder Stirling 
engine thermodynamic model is shown in figure 1. Each working space con- 
tains seven volumes - one each for the expansion space, the heater, the 
cooler, and the compression space, and three for the regenerator. A thermal 
time constant for the regenerator mesh is associated with each regenerator 
gas volume. Within each control volume, the continuity and energy equations 
are solved and a simplified first-order momentum term (flow resistance) is 
calculated. Upwind differencing is used to calculate interface volume tem- 
peratures (primed temperature variables in fig. 1) for use in the energy 
equation. All symbols shown in figure 1, as well as a complete symbols list, 
are given in appendix A. The equations used to model the engine are given 
in appendix B. 

A schematic of the drive dynamics is shown in figure 2. Differential 
forces on the pistons are translated into torque through the vehicle drive 
geometry and summed to form total indicated torque. Torque due to engine 
friction is subtracted to form brake torque. This torque is available to 
drive auxiliaries and the vehicle load. The vehicle inertia and gear ratio 
are used in computing the effective load. The net torque is integrated once 
to give engine speed (PSIDT) and once again to give crank angle (PSI). The 
crank angle is used to generate piston position (using the crank geometry). 
The model can be run in this manner, or piston position can be forced as a 
function of time with torque and cycle performance calculated at constant 
speed. 

A model of the working fluid supply system is shown in figure 3. This 
system is actually part of a mean pressure control system which modulates 
engine power by changing the amount of working fluid in the cycle. In order 
to accelerate the engine, working fluid is supplied to the engine. A slot 
in the piston rod is used to time the injection of working fluid into the 
various compression spaces in the engine. The timing is arranged to supply 
the fluid only when the pistons are near bottom (minimum) stroke. During 
this period, the pressure in the associated compression space is near a max- 
imum, and the injected fluid functions to increase rather than decrease the 
engine torque. The check valves (fig. 3) remain open as long as the supply 
pressure (PSOURC) is greater than the corresponding compression space pres- 
sure. Also included in the supply model are rod leakage effects. These are 
indicated by WSSCN in figure 3 and allow study of the effects of mistimed 
flow injection into the system. 
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USERS MANUAL 


Simulation Flow Diagram 


The overall simulation structure is shown in figure 4. Run conditions 
are set in MAINST. Subroutine ICSTIR is then called to calculate initial 
conditions. PISTIN is called from ICSTIR to calculate the initial piston 
position. ICSTIR then calls INTEGR, which is the integration subroutine. 
INTEGR handles the incrementing of time, data output, and run termination. 
INTEGR calls the Stirling engine simulation subroutine STIRSM to obtain the 
information it needs for the Jacobian matrix that is required by the 
backward-difference integration. STIRSM calls TORLOS to calculate engine 
friction, auxiliary losses, and vehicle load effects; PISTIN, to calculate 
piston positions from crank angle and crank geometry; and HYSUPY, to calcu- 
late flow into the engine, the phasing of the injected flow, and rod leakage 
effects. STIRSM then returns to INTEGR. 

In figure 4 the body of INTEGR is shown within a dashed line. The sub- 
routines enclosed are actually subroutines to INTEGR. While they are called 
by and return to INTEGR, they are shown as calling one another to illustrate 
the looping that takes place within the subroutine. Once a Jacobian matrix 
is calculated, INTEGR calls DMINV for a matrix inversion, or BRYDON for up- 
dates to a previously generated matrix. Subroutine BRYDON contains the 
Broyden update algorithm. This algorithm permits the simulation operating 
point to be changed without having to calculate a new Jacobian matrix. The 
algorithm does this by continually updating the known inverted Jacobian ma- 
trix. This can significantly reduce the computation time, since generating 
Jacobian matrices and inverting the matrices can be very time consuming. 

The use of the algorithm will be discussed in more detail later. 

Convergence is then checked. If the simulation is not converged, more 
passes are made through INTEGR to generate better iteration guesses or, if 
necessary, a new Jacobian matrix. If converged, INTEGR calls TRAP to do a 
trapezoidal integration on the product of the pressure and volume in both 
the expansion and compression spaces. TRAP returns to INTEGR, which then 
calls OUTSTR if a printout is desired at the current time. OUTSTR returns 
to INTEGR; then GUESSE is called to predict the next set of state variables 
for the next time step. This is done until the desired maximum number of 
time steps have been run, at which time INTEGR terminates the run. Flow 
charts for all subroutines are given in appendix C. 


Program Setup 


The program setup is done in the main program MAINST. Switches are set 
to indicate the type of transient desired. Engine geometry and flow data 
are input data. Tables I to XIII indicate the required input as well as the 
options available. In table I, baseline engine geometry data are listed. 

The total regenerator volume (VR) excludes the volume of the mesh. Tables 
II to IV list the required heater, cooler, and regenerator data, respective- 
ly. The initial guesses for the temperature distribution within a working 
space are shown in table V. All four working spaces start out with this 
distribution as a first guess. All of the simulation constants are listed 
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in table VI. Flow resistances between the volumes are shown in table VII. 
Hydrogen constants are listed in table VIII. Data for the implicit 
integration method are listed in table IX. In most cases, the indicated in- 
tegration settings should lead to convergence. However, experience has 
shown that increasing MPAS and/or T0LPC6 can be helpful for difficult con- 
vergence problems. 

Another method of overcoming convergence problems is to increase TOLSS. 
This is done internally in the program if the number of iterations using one 
matrix exceeds 20 passes, or if a new matrix needs to be generated. Once 
the program converges, TOLSS is set back to .0001 for the next time point. 

If MPAS is exceeded, the simulation will continue to run, but will produce a 
debug output at that time point indicating that MPAS has been exceeded. 

The different switch options available can be selected by a set of 
switches defined in table X. Table XI lists all the run conditions that 
must be specified, while table XII defines the load characteristics. The 
cycle data are defined in table XIII. Note that if a supply transient is 
not desired, NCYSUP and NCYSTP from table XI must be set greater than NUMBCY. 


Output Options 


The user may select from a number of output options and a debug option. 
If ICALC = 1, a short printout is specified. The short printout header is 
shown in table XIV. Output variables include the heater temperature (TWH), 
cooler temperature (TWO), cycle mean pressure (CYCLPR), engine drive mode 
(NONENG), number of cycles to be run (NUMBCY), piston rod leakage area 
(ALEAK), supply pressure (PSOURC), hydrogen supply temperature (TSOURC), the 
time point to start the hydrogen supply (ISUPST), and the time point to end 
the hydrogen supply (ISPSTP). 

Results are also printed at specified time steps. Power, torque, crank 
angle, mean pressure, and engine speed are printed. ITRAN indicates the 
number of time steps taken up to and including the printout; KWORK is a 
counter indicating the number of time steps in the calculation of the 
pressure-volume area (note that if KWORK = 201, 200 time steps were taken, 
since the first time step printed out is the initial condition). 

If ICALC = 0, a long printout is specified; the long printout header is 
shown in table XV. At each specified time step a complete listing of all 
engine variables is generated. Note that the output is arranged in rows of 
four (except for the overall engine parameters). The first row corresponds 
to the first working space, the second corresponds to the second, etc. The 
output variables include all pressures, temperatures, and stored masses for 
each working space, as well as the power, torque, and supply flow. The last 
row of the long printout contains the same overall engine parameters as the 
short printout. 

The user may select a debug option to aid in solving problems that may 
occur in the simulation. This option is specified in MAINST by NOBUG. If 
NOBUG = 1, the debug option will not be activated unless (a) there is a 
problem with convergence (the maximum number of iteration passes is exceeded 
without convergence), or (b) there is a problem in generating a partial de- 
rivative for the Jacobian matrix. The debug printout comes from STIRSM. 

The debug option is selected by the user by setting NOBUG = 0. A sample of 
the debug printout is shown in table XVI. This option is further described 
in appendix C. 
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OUTPUT - TEST CASE 


Supply Transient 


A supply transient is provided as an output test case. The assumed 
initial conditions for the transient are a mean pressure of 5 MPa and a 
speed of 2000 rpm. After 10 cycles, hydrogen is injected to accelerate the 
engine. The hydrogen supply pressure is 10 MPa and the supply temperature 
is 283 K. Piston rod leakage is assumed to be zero for this test case. The 
corresponding Fortran input for MAINST is shown in table XVII. The 
transient is defined by CYCLPR = 5, SPDRPM = 2000, NCYSUP = 10, NUMBCY = 100, 
PSOURC = 10, TSOURC = 10, and ALEAK = 0. The Broyden update algorithm is 
used (IBRYTH = 1), the step size is determined by the specified 200 points 
per cycle (NPTPCY = 200), and a short printout is desired every 200 points 
(ICALC = 1 and IPRTOP = 200). 

The printout for this test case is given in table XVIII. Note that 
time is incremented so that the specified number of integration points per 
cycle is obtained. Therefore, the time step is a function of engine speed. 
Sixty-one engine cycles (12 201 time steps at 200 per cycle) were simulated 
before the specified limit on CPU time was exceeded. On the IBM 3033, the 
simulation ran for 60 min of CPU time, taking about 1 min of CPU time per 
engine cycle. 

Figure 5 shows a comparison of acceleration results from this simula- 
tion (dashed lines) with results from the simplified model (solid lines). 

Data for power, torque, and mean pressure are plotted against time from the 
start of the hydrogen supply (both simulations were run for 10 cycles before 
supply input). At a 5 MPa mean pressure, the seven-volume model predicts 
more power (and thus more torque) than the simplified model (as seen by the 
steady-state maps of ref. 2). Thus, with the same load, losses, and gear 
ratios used on both simulations, the seven-volume model reaches steady-state 
at a higher speed than the simplified model. 

At time equal zero in figure 5, the seven-volume model is running at 
2400 rpm, while the simplified model is running at 2000 rpm. A hydrogen 
supply transient is then simulated in both cases, assuming a 10-MPa supply 
bottle. Pressure rises at about the same rate in both simulations. Power 
rises in both simulations at about the same rate, with the difference re- 
maining fairly constant. Total torque also rises in both simulations. 
However, the difference between the total torques decreases slightly because 
of the difference in speed-torque maps as determined by each of the models. 
Both simulations give basically the same dynamic response; however, the 
three-volume model uses much less CPU time (about 1 sec of CPU time per 
cycle) because of its fewer state variables (14 versus 70). 

Thus, both simulations are useful for controls analysis. The three- 
volume simulation (ref. 5) can be used for evaluating many control strate- 
gies or doing parametric studies, while the seven-volume simulation can be 
used to study detailed performance during specific controls transients. 
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Cycle Analysis 


The simulation can be used to look at individual cycles. Results from 
a steady-state case (no supply) are shown in figure 6. Piston position, in- 
dividual torques, total torque, compression space pressures, and compression 
space temperatures are shown. Speed is constant at 2000 rpm for this case. 
Results are similar to those presented for the simplified model in reference 
4, except for the overall torque and compression space temperatures. Over- 
all torque (fig. 6(c)), shows a low frequency oscillation at steady state 
which was not observed with the simplified (ref. 5) model. Further investi- 
gation, however, showed that the simulation had not quite reached steady 
state. Data are presented after eight engine cycles. The oscillations are 
approximately 1 N-m in 61 N-m, or about 1.6 percent. The reason for the 
torque oscillation is that the pressures within the working spaces are not 
at steady state, due to the slow response of the regenerator metal lumps 
(not included in the three-volume model). 

Figure 7 shows the supply transient response of the engine with no pis- 
ton rod leakage. The piston motion, individual torques, total torque, com- 
pression space pressures, compression space temperatures, supply flow, and 
rod leakage flow are plotted for three cycles. Again, results are similar 
to those presented for the simplified model in reference 4. Note, however, 
that compression space temperatures (fig. 7(e)) do drop because of the in- 
jection of cold working fluid and that compression space pressures (fig. 
7(d)) rise because of the increase in stored mass. No torque droop is indi- 
cated in figure 7(c). 

The same supply transient was run with a rod leakage area of 3.2 cm^ 
(0.5 in^); results are shows in figure 8. No torque droop is indicated in 
figure 8(c), even though there is a mistiming of the injected flow. Com- 
pression space temperatures, however, do fall more rapidly than in the non- 
leakage case. 


CONCLUSIONS 


A detailed seven-volume, four-working-space Stirling engine computer 
simulation has been developed. The simulation model includes drive dynamics 
and vehicle load effects. The simulation program includes subroutines which 
allow for simulation of control strategies, such as acceleration of the en- 
gine by adding hydrogen to the system. 

All input data required to run the program are described. Flow charts 
of the overall simulation and the subroutines are also given. The simula- 
tion is modular to allow for easy modification of the model. 

An input routine is provided in which the user specifies the transient 
to be run and supplies the necessary engine geometry data. The type of out- 
put is also selected by the user. Very detailed printouts at prescribed in- 
tervals can be selected; alternatively, less detailed printouts can be se- 
lected. The integration step size and the printout interval do not have to 
be the same. 

Simulation test cases have been run, and results have been compared 
with results obtained from a simplified three-volume simulation. The simu- 
lations produce similar results, with the more detailed simulation (with 
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temperature dynamics) exhibiting more engine output power than the simpli- 
fied simulation at the same engine speed and pressure. The detailed simula- 
tion shows variations in compression space temperature due to injected mass; 
the simplified simulation has no temperature variation. 

The detailed simulation can be used to predict steady-state or tran- 
sient performance of the engine. Engine accelerations can be run with or 
without piston rod leakage. The acceleration data presented herein can be 
used to verify the transient operation of the simulation program when actual 
engine data become available. 
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APPENDIX A 


NOMENCLATURE 


Engineering symbols 


A 

A 

A 

C 

C 

F 

G 

h 

I 

J 


d 

r 

P 

V 


i. 


P 


Q 

R 

S 

T 

t 


V 

Vel 

w 


2 

area, m'- 

p 

piston area, cm^- 
2 

rod area, cm*- 

specific heat at constant pressure, J/(kg-K) 
specific heat at constant volume, J/(kg-K) 
force, N 
gear ratio 

p 

heat transfer coefficient, J/(sec-m^-K) 
inertia, N-m-sec'^ 

mechanical equivalent of heat, 1.0 (N-m)/J 

piston stroke, m 

pressure, N/m^ 

heat flow rate, J/sec 

gas constant, (N-m)/(kg-K) 

connecting rod length, m 

fluid resistance, (N-sec)/{kg-m ) 

temperature, K 

time, sec 

torque, N-m 

•3 

volume, m*^ 

vehicle velocity, km/hr 
mass, kg 


* 1 . 1 ; * 1 , 2 ; 

'^1,3; *1,4 

* 2 , 1 ; * 2 , 2 ; 

*2,3; *2,4 


mass flow rate, kg/sec 

flow rate from expansion space to heater, kg/sec 
flow rate from heater to 1st regenerator, kg/sec 
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'^3,1; *3,2; 
*3,3; *3,4 

*4,1; *4,2; 
*4,3; *4,4 

*5,1; *5,2; 
*5,3; *5,4 

* 6 , 1 ; * 6 , 2 ; 
*6,3; *6,4 

Y 

L 


Subscripts; 

a 

c 

d 

e 

f 

h 

i 

in 

J, k 
tn 

out 

P 

r 

rr 

s 

sup 

vol 

w 

0 


flow rate from 1st regenerator to 2nd regenerator, kg/sec 
flow rate from 2nd regenerator to 3rd regenerator, kg/sec 
flow rate from 3rd regenerator to cooler, kg/sec 


flow rate from cooler to compression space, kg/sec 

ratio of specific heats 
change in 
angular position 

auxiliaries 

cooler 

drag 

engine 

mechanical friction 
heater wall 

position in working space 
into 

working space 

mesh 

out of 

piston 

rod 

rolling resistance 

stored 

supply 

gas volume 

wheel 

dead volume 


Superscripts; 

I 


intervolume 

derivative 
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Computer Variables 


AD 

piston area, cm^ 

ADSl, ADS2,) 



slot area in piston rods, cm^ 

ADS3, ADS4 J 


AJ 

mechanical equivalent of heat, 1.0 (N-m)/J 

ALEAK 

piston rod leakage area, cm^ 

ALPHA 

crank angle lag, deg 

AR 

2 

piston rod area, cm 

AWC 

cooler heat transfer area, cm^ 

AWH 

2 

heater heat transfer area, cm 

AWR 

2 

regenerator heat transfer area, cm 

BIGNUM 

scalar for matrix predictor (biggest number) 

CMCT 

scratch matrix in BROYDEN algorithm 

CP 

specific heat at constant pressure, J/(kg-K) 

CVM 

specific heat of the mesh, J/(kg-K) 

CYCLPM 

maximum cycle pressure, MPa 

CYCLPR 

current cycle pressure, MPa 

DEGR 

conversion factor, rad/deg 

DELT 

time change, sec 

DELTAV 

vector change in guess variables 

TO( 

vector change in guess variables (BROYDEN algorithm) 

DELY 

vector change in error variables (BROYDEN algorithm) 

E 

error vector 

EMAt 

jacobian matrix for 70th order system 

mmi 

past error vector 

FDRAG 

force due to drag, N 

FORCl, F0RC2,) 



piston forces, N 

F0RC3, F0RC4 j 


FRAC 

external control for matrix convergence 

FRRF 

force due to rolling resistance, N 

G 

O 

gravitational constant, kg/(MPa-sec ) 

GAMMA 

ratio of specific heats 
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GR 

GRAT 

GTRAN 

HC 

HH 

HR 

IBRYTH 

ICALC 

ICON 

IHPCNV 

INDICA 

I PRIOR 

I ROT 

ISPSTP 

ISS 

ISUPST 

I TRAN 

ITRMAX 

KBROY 

KWORK 

KWRIT 

MATRIX 

MATTOT 

MPAS 

N 

NCYSTP 

NCYSUP 

NITER 

NMAT 

NOBUG 

NONENG 


gear ratio 
overall gear ratio 
transmission gear ratio 

2 

cooler heat transfer coefficient, J/(sec-cm -K) 

2 

heater heat transfer coefficient, J/(sec-cm -K) 
regenerator heat transfer coefficient, 

J/(sec-cm^-K) 

switch for BROYDEN algorithm 

switch for output printout 

counter for converged errors 

switch for matrix generation at every point 

switch for sign change on guess variable 

switch for the number of printouts per cycle 

switch for indicating a complete crank rotation 

time point for stopping supply 

switch for initial conditions 

time point for start of supply transient 

counter for time steps 

maximum number of stime steps 

counter for calls to the BROYDEN algorithm 

number of points in pressure-volume integration 

counter for printout 

switch for generating a new Jacobian matrix 
counter for the number of matrices generated during 
a transient 

maximum allowable iteration passes 
system order 

cycle at which to end supply 
cycle at which to start supply 
counter for the number of iterations at a time point 
counter for the number of Jacobian matrices gene- 
rated at a time point 
switch for a debug printout 
type of piston motion 

0 forced 

1 calculated 
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NPTPCC 

NPTPCY 

NSTP 

NTMAX 
NUMBCY 
PAUX 
PAVEMP 
PCI, PC2,| 

PCS, PC4 I 

PCNCHG 

PCOLDl, PCOLD2, 
PCOLD3, PC0LD4 , 
PDRAG 
PENG 

PEI, PE2,] 

PE3, PE4 / 
PFRICT 


number of integrations per cycle plus one 

number of integrations per cycle 

switch for storing past converged scale factors for 

iteration guesses 

maximum number of iteration variables 
number of engine cycles to run 
power loss due to auxiliaries, MPa 
cycle pressure, MPa 

gas pressure in compression spaces, MPa 

iteration convergence rate 

gas pressure in cooler spaces, MPa 

power loss due to vehicle air drag, kW 

total power loss due to friction and auxiliaries, kW 

gas pressure in expansion spaces, MPa 

power loss due to engine friction, kW 


PHI, PH2,j 

PH3, PH4 j 

PIE 

PLOAD 

PNET 

POSDEG 

POSDEO 

POWER! 

PRRF 

PRll, PR12,' 
PR13, PR14 J 
PR21, PR22, 
PR23, PR24 J 


gas pressure in the heater spaces, MPa 

constant (3.14167) 

total power required by the load, kW 

net power, kW 

crank angle, deg 

initial condition crank angle 

total engine power output, kW 

power loss due to rolling resistance, kW 

gas pressure in first regenerator spaces, MPa 
gas pressure in second regenerator spaces, MPa 
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PR31, PR32,| 
PR33, PR34 ) 

PSI 

PSIDDT 
PS IDT 
PSOURC 
PWRl, PWR2,| 

PWR3, PWR4 / 

QCl, QC2,| 

QC3, QC4j 

QHl, QH2.| 

QH3, QH4 I 

R 

RANK 

RAT 

RATIO 

REF 

RODL 

RST 

RWHEEL 

SES 

SESP 

SPDMAX 

SPDRPM 

STROKE 

TAUX 

TCO 

TCI, TC2,. 

TC3, TC4 f 
TCOLDO 

TCOLDl, TCOLD2, 
TCOLD3, TC0LD4 


gas pressure in third regenerator spaces, MPa 
crank angle, rad 

2 

crank angle acceleration, rad/sec 
crank angle velocity, rad/sec 
hydrogen supply pressure, MPa 

power out of cylinders, kW 


cooler heat flow rate, J/sec 


heater heat flow rate, J/sec 

gas constant, (MPa-cm )/(K-kg) 

change centigrade to kelvin 

scaler change on step size for iteration 

largest step size change 

desired value of the summation of errors 

crank rod length, cm 

flow resistance, (MPa-sec) /kg 

radius of car wheel, cm 

summation of squares of present errors 

summation of squares of past errors 

maximum engine speed, rpm 

current engine speed, rpm 

piston stroke length, cm 

torque loss due to engine auxiliaries, N-m 

initial compression space temperature, K 

gas temperature in compression spaces, K 
initial cooler space temperature, K 
gas temperatures in cooler spaces, K 
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TDRAG 

TEMP 

TEMPI 

TEMP2 

TTRP3 

TENG 

TEO 

TEl, TE2,| 

TE3, TE4 ) 

TFRICT 

THO 

THl, TH2, 
TH3, TH4 


torque loss due to vehicle air drag, N-m 

Broyden update scalar 

Broyden update scalar 

Broyden update vector 

Broyden update vector 

total torque loss due to engine auxiliaries and 
friction, N-m 

initial expansion space temperature, K 
gas temperatures in expansion spaces, K 

torque loss due to engine friction, N-m 
initial heater temperature, K 

gas temperature in heater spaces, K 


TIME 

TLOAD 

TNET 

TOLPCG 

TOLSS 

TOLl 

T0L2 

TORQT 

TORQl, T0RQ2, 
T0RQ3, T0RQ4 


time, sec 

total torque loss required by the load, N-m 
net torque, N-m 

convergence rate at which a decision is made to 
generate a new Jacobian matrix 
error tolerance 

lower limit for good partial derivative 
upper limit for good partial derivative 
total engine torque, N-m 

individual torques, N-m 


TRRF 

TRIO 

TRll, TR12, 
TR13, TR14 J 


torque loss due to rolling friction, N-m 
initial regenerator gas temperature, K 

first regenerator space temperature, K 


TR20 

TR21, TR22, 
TR23, TR24 
TR30 


initial regenerator gas temperature, K 
second regenerator space temperature, K 
initial regenerator gas temperature, K 


15 



TR31, TR32,| 
TR33, TR34 ) 
TRQl, TRQ2,| 
TRQ3, TRQ4 ) 


third regenerator space temperature, K 


individual piston torques, N-m 


TSOURC 

TWC 

TWH 

TWRll, TWR12 
TWR13, TWR14 
TWR21, TWR22,| 
TWR23, TWR24 J 
TWR31, TWR32,| 
TWR33, TWR34 J 

Tl,l; Ti.2; 

► 

Ti,3; Ti,4 , 
T2,1; T^.2; 

► 

T2,3; T2.4 , 


working fluid source temperature, K 
cooler wall temperature, K 
heater wall temperature, K 

mesh temperature in first regenerator segments, K 
mesh temperature in second regenerator segments, K 


mesh temperature in third regenerator segments, K 


intervolume temperature between expansion space and 
heater space, K 


intervolume temperature between heater space and 
first regenerator space, K 


^3,1; T3,2;' 
■^3,3; T3,4 , 

Ti,l; n, 2 ; 

T4,3; T4,4 
T5,1; T5,2; 
■^5,3; T5,4 , 


intervolume temperature between first regenerator 
space and second regenerator space, K 


intervolume temperature between second regenerator 
space and third regenerator space, K 


intervolume temperature between third regenerator 
space and cooler space, K 
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intervolume temperature between cooler space and 
compression space, K 


VCl, VC2, 

VC3, VC4 
VCOLD 

vcMv 

VDELTA 

VDDt 

VDOTT 
VEl, VE2,j 

VE3, VE 4 J 

VGUESS 

VH 

VKPH 

VMST 

VOC 

VOE 

VR 

VS 

VSM 

W? 

WSCl, WSC2,| 
WSC3, WSC4 / 
WSCLDl, WSCLD2, 
WSCLD3, WSCLD4J 
WSEl, WSE2,| 
WSE3, WSE4 / 
WSHl, WSH2,) 
WSH3, WSH 4 / 


compression space volumes, cm^ 
cooler volume, cm^ 

vector of converged state variables from previous 
time step 

initial perterbation of state variables 

vector of state variable derivatives at current time 

step 

vector of average value of state variable derivatives 
expansion space volumes, cm^ 

initial guess vector at each time step 
3 

heater volume, cm 
car velocity, km/hr 

vector of changes in guess variables during iteration 

3 

compression space dead volume, cm 

3 

expansion space dead volume, cm 

3 

regenerator volume minus mesh, cm 
guess vector during iteration 
vector of saved conveyed guess variables 
state vector 

compression space stored masses, kg 


cooler stored masses, kg 


expansion space stored masses, kg 


heater space stored masses, kg 
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WSM 

WSOl, WS02,| 

WS03, WS04 I 

WSRll, WSR12,| 

WSR13, WSR14 I 

WSR21, WSR22,| 

WSR23, WSR24 1 

WSR31, WSR32A 

WSR33, WSR34 ) 

WSSl, WSS2, 

WSS3, WSS4 

WSSCNl, WSSCN2,| 

WSSCN3, WSSCN4 J 

WSTOTl, WSTOT2,] 

WST0T3, WST0T4 I 

WTl, WT2,| 

WT3, WT4 I 

WTCAR 
WTENG 
WTOT 
WTWHEL 
XDl, XD2,| 

XD3, XD4 I 

XD(3) 

XD30 

XXX 

YW 

im 


mass of the mesh, kg 
intermediate supply flows, kg 


first regenerator stored masses, kg 


second regenerator stored masses, kg 


third regenerator stored masses, kg 


supply flows, kg/sec 


rod leakage flows, kg/sec 


total supply flows, kg 


stored mass in each working space, kg 

mass of the car, kg 
mass of the engine, kg 
total stored mass, kg 
mass of the wheels, kg 

piston positions, cm 

piston 3 position, cm 

piston 3 starting position, cm 

summation of squares of the changes in the errors to 

the max error 

scalar vector for changes in guesses 
Jacobian matrix for 68th order system 
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APPENDIX B 


MODEL EQUATIONS 


The equations used to model the fluid dynamics and thermodynamics are 
ordinary differential equation approximations of the complex partial deriva- 
tive flow equations. Mass flow is assumed to occur due to pressure differ- 
ential between the gas nodes. In figure 1 the pistons are numbered in the 
order in which they reach top stroke, as is indicated by the arrows in the 
piston heads. A double subscript notation has been adopted in which the 
first subscript denotes a volume position within a working space, and the 
second subscript denotes the working space. 

Representative equations are 



(P 


^i,j»i+l,j 


) 


j = 1,2. ..4 
i = 1,2... 6 


where the arrow indicates flow resistance measured between the volumes. 




and 


P 


i.j 




j = 1,2. ..4 
i = 0,1. ..6 




0.0 


j = 1,2... 4 
i = 1,2. ..7 


Volumes are assumed to vary with piston position. For volume 1 in working 
space 1 



Piston position is a function of crank angle: 
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The energy equation used to calculate the change in temperature in a volume is 
complicated by the oscillatory nature of the Stirling cycle. The form of the 
equation is 


w 


Q . Work. . 

T' = ft - fvT' , - T. .1 - ft (yT! - T .) + ' ^ 

l,j 1,0^ i.j 


where 


j = 1,2.. .4 
i = 1,2.. .7 


and where the heat flows into and out of the gas are modeled as 


Q = hA(AT) 


Primed temperature variables denote interface volume temperatures and are 
determined by upwind differencing (ref. 1). The method assumes that the gas 
temperature has the same temperature profile as the regenerator matrix tem- 
perature. Representative equations for positive and negative flows are 


ft. . > 0, T! . = T. . 
i,J - ’ i.J i.J 


) w,.At(T^ - T„ ) 
"^3,j ^ \ *"3,j *"5,j 

47^ 4.0 w. 


ft, , < 0, T! , = T 




i.j “ 1+1, j 


(Tm - ) 

"^3,3 ^ ^ "*3,j 


"O" 


4.0 w. 


M+l,j 


j = 1,2. ..4 
2 < i < 5 


ft. . > 0, T'. . = T. . 
I.J- i.J i.J 


ft. . < 0, T'. . = T. . 
I.J I.J 1+1. J 


i = 1,6 


The regenerators are modeled as thermal lags; 


h.A. 

T = (T. . - T ) 

"’i i S '^s "’i i 

I.J Pm s i,j 


j = 1,2. ..4 
i = 3,4,5 


Torque is calculated as a function of differential forces on the pis- 
tons, which are summed through the drive geometry. A representative equa- 
tion is 
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Also included in the model are simplified vehicle load effects and en- 
gine power losses due to mechanical friction and auxiliaries. Figure 2 
shows a schematic of how these losses are calculated. Torques from the pis- 
tons are summed to form indicated torque. Torque due to mechanical friction 
(TFRICT) is subtracted to form brake torque. The resultant torque is avail- 
able for engine auxiliaries (TAUX) and also for vehicle load effects, such 
as rolling resistance (TRRF) and drag (TDRAG). 
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APPENDIX C 


FLOW CHARTS 


This appendix contains flow charts for the main program and all the sub- 
routines of the simulation. Most flow charts are straightforward, with the 
possible exception of that for subroutine INTEGR. INTEGR performs both the 
incrementing of time and the integration of the state equations. The inte- 
gration scheme is implicit» It is a backward-difference integration which 
uses a multivariable Newton-Raphson iteration for convergence at a time 
point. This type of integration is stable for both large and small step 
sizes. This stability is important when there is a large spread in eigen- 
values for the system, which is the case when both pressure-flow and temper- 
ature dynamics are simulated. The implicit integration scheme can handle 
the widespread dynamics while insuring stability. To help one understand 
how subroutine INTEGR works, the following description is provided. To help 
one follow the INTEGR flow chart, statement numbers corresponding to the 
Fortran listing are given in the flow chart. 


Integration Scheme 


The integration scheme is a backward difference method which uses a mul- 
tivariable Newton-Raphson iteration scheme for convergence. Guess variables 
are updated by using the old guess vector V?old* the current error vector 
r, and an inverted Jacobian matrix EMAT : 


new 


-1 

- EMAT X E + VS 

old 


(Cl) 


EMAT is a Jacobian matrix of partial derivatives (i.e., changes in error 
variables with respect to changes in state variables): 


EMAT(I,J) = [E(J) - ERRBSE(J)]/DELTAV(I) (C2) 


Updating takes place when the errors are converged within tolerance. With 
this technique, both steady-state and transient solutions can be obtained by 
redefining the error variables. In steady state all states are at rest; thus 

VDDT=0.0 = E' (C3) 
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For a transient case 


VDOT X DELT - = E (C4) 

Equations (C3) and (C4) are converged when all the elements of F are with 
a specified tolerance. It should be noted that all the equations are scaled 
by the converged state to help convergence. 

In order to use this method, a Jacobian matrix must be calculated 
(usually by finite differences) and then inverted. This can be very time 
consuming. Therefore, the logic in INTEGR allows for calculation of a new 
matrix only under adverse conditions. Examples might be (1) too low a con- 
vergence rate of the simulation (PCNCHG < TOLPCG), or (2) the number of 
allowable passes (MPAS) being exceeded. Also to speed convergence, the error 
tolerance (for eq. (C4)) is increased if more than one matrix must be calcu- 
lated at the same time or if more than 20 passes are required using the same 
matrix. Once convergence is obtained, the tolerance is reset to its nominal 
value. 


Perturbation Calculation 


Several features in INTEGR help the implicit integration scheme converge. 
Of primary importance is the generation of a "good" Jacobian matrix. All 
the partial derivatives must be representative of the linear behavior of the 
system at a given operating point. Since finite differences are used, the 
sizes of the perturbations of the states are important: If the perturba- 

tions are too large, errors will be introduced by the system nonlinearities; 
if they are too small, the partial derivatives will be in error due to 
numerical problems (without double precision arithmetic). 

Thus, a tuning mechanism has been introduced into INTEGR to optimize the 
sizes of the perturbation. First the sum of squares of all the changes in 
the errors is calculated for each perturbation. Once this is done, the 
"goodness" of the partial is checked by calculating 


XXX 


1 

N 



[E(I) - ERRBSE(I)]^ 


for each state variable and then checking if 


TOLl < XXX < T0L2 


(C5) 

(C6) 


If all XXX' s fall within the tolerance band, the matrix is considered good. 
For this simulation, TOLl = 0.001 and T0L2 = 0.01. Since all the errors 
are scaled, this tolerance band lies between 0.1 and 1 percent. For a more 
linear system, the band could be larger; for a more nonlinear system, 
smaller. 
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Scaling of Perturbations 


In general, for the initial perturbations at a point, the XXX's will 
not fall in the tolerance band described above. Thus, INTEGR scales the 
perturbations to try to force the XXX's within the band. This is done by 
calculating 

YYY = REF/XXX (C7) 

for each state variable. REF is defined as being the center of the tolerance 
band: 


REF . 1T0L1^^^T0L2) ,^8 

Once a set of YYY's has been calculated such that the XXX's fall within 
the band, the set of YYY's is stored. After this has been done for all N 
states, the scaling vec tor YYY is generated. When a new matrix is needed, 
the scaling vector YYY is applied to the current states to determine first 
guesses for the perturbations needed to obtain new partial derivatives. If 
for any state variable the new XXX falls outside the tolerance band, YYY 
is updated and the new result stored. This method generally reduces the 
number of passes required for subsequent matrix generation. 


Error Messages 


In generating a partial derivative, a situation may arise where XXX 
never gets within tolerance. When this happens, the program prints an error 
message: "CHECK INPUT - BAD PARTIAL DERIVATIVE", prints a debug output to 
help the user diagnose the problem, and then stops the simulation. This is 
the only time when the simulation is stopped, except for the normal exit 
(i.e., ITRAN incremented to its maximum value, ITRMAX). In general, noncon- 
vergence will occur when inconsistent coding is added to the simulation. 

One example would be if a piston ring leakage flow was calculated and sub- 
tracted from one working space but not added to the adjacent working space. 

Another error message occurs when the simulation does not converge. 

This situation occurs when MPAS (set at 50) is exceeded. A message is 
printed; for example, "ITERATION FAILURE 69 51 32 15". The numbers printed 
out are the number of converged errors (may be any number from 0 to N“l; 

69 is shown here), the number of iteration passes (MPAS + 1), the number of 
times the BROYDEN update algorithm has been called for the current matrix 
(32 here), and the point at which the convergence failed (ITRAN). In this 
situation, a debug output is printed from STIRSM which indicates 

I counter for the system order 

VS current guess variables 

VCONV past converged guess variables 

VWS state variables 

VDOT current state derivative 

VDOTT converged state derivative between current time point and past con- 
verged value 
E current errors 
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After the printout, the simulation continues. Note that with this method a 
convergence failure may occur, but the errors will be very close to the 
tolerance band. After the failure, the simulation may recover. For this 
reason the simulation is allowed to continue and the user may make a judg- 
ment as to the validity of the data after a convergence failure. 

The occurrence of many convergence failures in a transient, however, 
usu- ally indicates a need for the user to increase tolerance or check his 
input and coding. 


Broyden Update Algorithm 


As stated earlier, it is time consuming to calculate a Jacobian matrix, 
and even more time consuming to invert a matrix for a large system. Usually 
double precision is required (as in DMINV). Lack of speed of the matrix in- 
version algorithm can be prohibitive, especially in a controls model where 
long transients are needed to evaluate controls schemes. One possible means 
of avoiding this problem is to generate the inverse, and then to continually 
update it with information gained as the inverse is used. One way to do 
this is by using the Broyden update algorithm, found in subroutine BRYDON. 
The use of the algorithm is controlled by IBRYTH, which is set in MAINST. 
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Main program MAINST 


Specify engine geometry and J 
heater, cooler, and 
regenerator heat transfer 
characteristics 


Specify flow resistances, 
hydrogen gas constants, 
matrix convergence dat^ 
temperature distribution, 
and constants 


Specify switch settings, run ; 
conditions, load conditions, 
and cycle data 


Calculate the number 
of time steps, the start and 
stop of the supply.and the 
InertlaSi set regenerator 
mesh temperatures 
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Subroutine BRYOON 


Subroutine OMINV 




Subroutine GUESSE 
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Subroutine IMEGR 
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Subroutine INTEGR - Continued. 


© © 
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Subroutine INTEGR - Continued, 



(D © 



© © 


© ® 
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Subroutine IMIEGR - Continued. 
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SubroutiM INTCGR - Continued. 
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Subroutine INTEGR • Continued. 
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Subroutine INTEGR - Continued. 
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Subroutine INTEGR - Concluded. 



1500 ) 
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Subroutine STIRSM. - Continued, 



Subroutine TORLOS 



Subroutine TRAP 
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TABLE I. - PROGRAM INPUT: ENGINE GEOMETRY 


Name 

Setting 

Description 

AD 

AE 

RODL 

VR« 

VOE 

VOC 

VH 

VCOLD 

STROKE 

ALPHA 

23.7613 

I. 1290 
10.0 

115.332 

II. 459 
32.5939 
52.4058 
29.2378 

4.0 

90.0 

Piston area 

Piston rod area 

Crank rod length 

Regeneration volume 

Expansion space dead volume 

Compression space dead volume 

Heater volume 

Cooler volume 

Piston strike 

Crank angle lag 


^Excludes mesh volume. 


TABLE II. - PROGRAM INPUT: HEATER DATA 


Name 

Setting 

Description 

HH 

71.4061 

Heater heat transfer coefficient 

AUH 

426.9669 

Heater heat transfer area 

THH 

705.0 

Heater wall temperature 


TABLE III. - PROGRAM INPUT: COOLER DATA 


Name 

Setting 

Description 

HC 

51.8739 

Cooler heat transfer coefficient 

AWC 

804.1274 

Cooler heat transfer area 

TWC 

86.0 

Cooler wall temperature 


TABLE IV. - PROGRAM SETUP: REGENERATOR DATA 


Name 

Setting 

Description 

HR 

168.9762 

Regenerator heat transfer coefficient 

AUR 

66870.834 

Regenerator heat transfer area 

CVM 

142.1345 

Specific heat of the mesh 

USM 

0.6523 

Mass of the mesh 
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TABLE V. - PROGRAM SETUP: GAS TEMPERATURE DISTRIBUTION GUESS 


Name 

Setting 

Description 

TEO 

624.44 

Expansion space temperature 

THO 

632.78 

Heater temperature 

TRIO 

505.56 

1st regenerator segment temperature 

TR20 

360.0 

2nd regenerator segment temperature 

TR30 

214.44 

3rd regenerator segment temperature 

TCOLOO 

105.0 

Cooler temperature 

TCO 

105.0 

Compression space temperature 


TABLE VI. - PROGRAM SETUP: CONSTANTS 


Name 

Setting 

Function 

RANK 

273.0 

Converts centigrade to Kelvin 

DEGR 

57.296 

Converts degree to radians 

PIE 

3.1416 

R 

4125.6 

Gas constant 

G 

10017.0 

Gravitational constant 

AJ 

1.0 

Mechanical equivalent of heat 


TABLE VII. - PROGRAM INPUT; FLOW RESISTANCES BETWEEN VOLUMES 


Name 

Setting 

Description 

RST(l) 

0.38010 

Resistance between expansion space 



and heater 

RST(2) 

0. 76020 

resistance between heater and 1st 

RST(3) 

0.76020 

regenerator segments 

resistance between 1st regenerator 

RST{4) 


and 2nd regenerator segments 

0.76020 

resistance between 2nd regenerator 

RST(5) 


and 3rd regenerator segments 

0. 76020 

resistance between 3rd regenerator 

RST(6) 


segment and cooler 

0.38010 

resistance between cooler and 



compression space 
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TABLE VIII. - PROGRAM INPUT: HYDROGEN CONSTANTS 


Name 

Setting 

Description 

CP 

GAMMA 

1292.13 

1.39 

Gas constant 
Gas constant 


TABLE IX. - PROGRAM INPUT: MATRIX INPUT DATA 


Name 

Setting 

Function 

VDELTA 

0.01 

Perturb Initial guesses by 1 percent 

FRAC 

1.0 

External control of iteration step 
magnitude 

TOLl 

0.001 

Lower limit on tolerance error for 
matrix linearity 

T0L2 

0.01 

Upper limit on tolerance error for 
matrix linearity 

TOLSS* 

0.0001 

Solution tolerance 

N 

68 

System order NONENG » 0 


70 

System order NONENG • 1 

NTMAX 

70 

Largest system order 

MPAS 

so 

Limits convergence failure to SO 
passes 

TOLPCd’ 

0.1 

Recalculate Jacobian matrix when 
convergence rate (PCNCHG) falls 
below TOLPCG - no BROYDEN 
algorithm. 


®If the number of Iterations Is greater than 20 or the number of 
matrices generated Is greater than 1, TOLSS Is set to 0.001 at 
that point. 

*>If BROYDEN is used, TOLPCG - 0. 
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TABLE X. ^ PROGRAM INPUTS: SWITCHS 


Name 

Setting 

Function 

ISS 

0 

Set up initial conditions 


1 

Run transient (set internally in 
program) 

I CALC 

0 

Detailed printout 


1 

Shortened printout 

MATRIX 

1 

If initial guesses do not converge 
the simulation, generates a 
Jacobian matrix 

NONENG 

0 

Force crank angel as a function of 
time 


1 

Run the simulation as an engine 

I PR TOP 

NN 

Printout data every NN points 

IBRYTH 

0 

Oo not use BROYDEN update algorithm 


1 

Use BROYDEN update algorithm 

IHPCNV 

0 

Use logic in program to determine 
need for a new Jacobian matrix 


1 

Generate a new Jacobian at every 
time point 

NOBUG 

0 

No debut output 


1 

Internally set to give a debug out- 
put when a convergence failure or 
matrix generation problem occurs 


TABLE XI. - PROGRAM INPUT: RUN CONDITIONS 


Name 

Setting 

Function 

ALEAK 

0 

Piston rod leakage area 


15.0 

Hydrogen bottle source pressure 

TSOURC 


Hydrogen bottle source temperature 



Desired engine speed 

CYCLPR 


Initial engine pressure 

SPDMAX 


Maximum engine speed 

CYCLPM 

15.0 

Maximum engine pressure 

NCYSUP* 

10 

Start hydrogen supply 

NCYSTP* 

200 

End hydrogen supply 


3If no supply transient is desired, make NCYSUP and NCYSTP 
greater than NUMBCY. 
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TABLE XII. - PROGRAM INPUT: LOAD CONDITIONS 


Name 

Setting 

Function 

GTRAN 

GR 

WTEN6 

UTWHEL 

WTCAR 

RWHEEL 

1.0/2.53 
1.5 
0 kg 

29.48 kg 
1420 kg 

30.48 cm 

Transmission gear ratio 
Gear ratio 
Engine mass 
Mass of a wheel 
Mass of a car 
Tire wheel radius 


TABLE XIII. - PROGRAM INPUT: CYCLE DATA 


Name 

Setting 

Function 

NPTPCY 

200 

Number of integration points/cycle 

NUMBCY 

3 

Number of desired cycles 
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TABLE XIV. - SHORT PRINTOUT (ICALC = 1) 

STIRLING ENGINE FOUR CYLINDERr SEVEN VOLUHES PER CVl.INDERf CONTROLS MODEL 
RUN CONDITIONS FDR THIS TRANSIENT 


TUH = 

978.3 

TWC = 

359.3 

CYCLPR = 

5.000 

NONENG 

0 

NUMBCY ^ 

100 

ALEAK 

= 0.0000 

PSOURC 

- 10.00 . 

TSOURC - 

283.3 

ISUPST = 

2001 

ISPSTP = 

40001 


TIME 

POUERT 

PFRICT 

PAUX 

PRRF 

PDRAG 

PLOAD 

PENG 

PNET 

ITRAN 

PAVEMP 

POSDEG 

TOROT 

TFRICT 

TAUX 

TRRF 

TDRAG 

TLOAD 

TENG 

TNET 

KUORK 

SPDRPM 


VKPH 



TABLE XV. - LONG PRINTOUT (ICALC = 0) 


STIRLING ENGINE FOUR CYLINDER. SEVEN VOLUMES PER CYLINDER. CONTROLS MODEL 


RUN 
TUH = 

CONDITIONS 

978.3 

FOR THIS TRANSIENT 
TUC = 359.3 

CYCLPR = 

5.000 

NONENO 


0 

NUMBCY = 

100 

ALEAK - 

0.0000 

PSOIJRC 

^ 10.00 

TSOURC - 

283.3 

ISUPST 


2001 

ISPSTP = 

40001 

TIME 

XDl 

XD2 

XD3 

XD4 

POSDEG 

SPDRPM 

PAVEMP 

UTOT 

DELT 

ITRAN 

USSl 

PEI 

PHI 

PRll 

PR21 

PR31 

PCOLDl 

PCI 

UTl 

PURI 

TRQl X 

USSCNl 

TEl 

THl 

TRll 

TR21 

TR31 

TCOLDl 

TCI 

THRU 

TUR21 

TUR31 

WSS2 

PE2 

PH2 

PR12 

PR22 

PR32 

PC0LD2 

PC2 

UT2 

PUR2 

TRQ2 

USSCN2 

TE2 

TH2 

TR12 

TR22 

TR32 

TC0LD2 

TC2 

TUR12 

TUR22 

TWR32 

USS3 

PE3 

PH3 

PR13 

PR23 

PR33 

PC0LD3 

PC3 

UT3 

PUR3 

TRQ3 

WSSCN3 

TE3 

TH3 

TR13 

TR23 

TR33 

TC0LD3 

TC3 

TUR13 

TUR23 

TUR33 

USS4 

PE4 

PH4 

PR 14 

PR24 

PR34 

PC0LD4 

PC4 

UT4 

PUR4 

TRQ4 

WSSCN4 

TE4 

TH4 

. TR14 

TR24 

TR34 

TC0LD4 

TC4 

TUR14 

TUR24 

TUR34 

TIME 

POWERT 

PFRICT 

PAUX 

PRRF 

PORAG 

PLOAD 

PENG 

PNET 

ITRAN 

PAVEMP 

PDSDEG 

TORQT 

TFRICT 

TAUX 

TRRF 

TDRAG 

TLOAD 

TENG 

TNET 

KUORK 

SPDRPM VKPH 



TABLE XVI. - DEBUG PRINTOUT (FROM STIRSM) 


1 VS VCONV VUS VDOT VDOTT E 

1 610.25000 0.15300290E-03 0 . 15981176E-03-0 .79999983E-01-0 . 39999992E-0I-0 .83716273E-0X 

2 608.25000 0 . 12<t55953E-03 0 . 12967623E-03-0 . 199999S1E-01-0 . 99999905E-02-0 . 53120777E-01 

3 603.25000 0 . 1 062997 9E-03 0 . 10975669E-03 0 . 19999981E-01 0 . 99999905E-02-0 . 1890 9215E-01 

9 599.25000 0 . 13073008E-03 0 . 13908699E-03 0 . 19999985E-01 0 . 99999905E-02-0 . 19199920 E-01 

5 596.25000 0 . 16979067E-03 0 . 17322700 E-03 0.00000000 0.00000000 -0 .20539172E-01 

6 593.25000 0 . 16693657E-03 0 . 16 900099E-03 0 .59999999E-01 0 . 29999997E-01 0.11632S19E-01 

7 593.25000 0 . 70013598E-03 0 . 710 92020E-03 0.00000000 0.00000000 -0 . 15903751E-01 

8 1616.0000 1616.0000 1616.0000 , -315990.56 -157795.25 -0 . 19692199E-01 

9 1631.0000 1631.0000 1631.0000 169707.06 82353.500 0 .75738952E-02 

10 1902.0000 1902.0000 1902.0000 358807.81 179903.88 0 . 19199929E-01 

11 1190.0000 1190:0000 1190.0000 259690.13 129820.06 0 . 17081585E-01 

12 878.00000 878.00000 878.00000 126617.88 63308.938 0 . 10815877 E-01 

13 681.00000 681.00000 681.00000 87375.000 93687.500 0 . 96227 979E-02 

19 681.00000 681.00000 681.00000 0.00000000 0.00000000 0.00000000 

15 1902.0000 1902.0000 1902.0000 0.00800000 0.00000000 0.00000000 

16 1190.0000 1190.0000 1190.0000 0.00800000 0.00000000 0.00000000 

17 878.00000 878.00000 878.00000 0.00000000 0.00000000 0.00000000 

18 650.55059 0 . 28928765E-03 0 . 28997737E-03-0 . 39999999E-01-0 . 20000 OOOE-01 0 .62577166E-02 

19 699.55059 0 . 19103967E-03 0 . 10972896 E-03 0.88000800 0.00000000 0.22199923 

20 697.55059 0 . 12036909E-03 0 . 11781685E-03-0 .20000000E-01-0 . 99999979E-02 0 .87002991E-02 

21 699.55059 0 . 19802665E-03 0 . 19922277E-03-0 . 39999966 E-01-0 . 1 9999981E-01 0 .59306997E-02 

22 639.55059 0 . 1921986 1 E-03 0 . 18580700 E-03-0 . 90000021E-0 1-0 . 20000011E-01 0 . 1769636 lE-0 1 

23 632.55059 0 . 18895737E-03 0 . 18019609E-03-0 . 99999969E-01-0 .99999982E-01 0 . 90395558E-02 

29 626.55059 0 .97153098E-03 0 .99658361E-03 0.23999995 0.11999995 0 . 91079712E-01 

25 1616.0000 1616.0000 1616.0000 -88617.063 -99308.531 -0 . 91127950 E-02 

26 2058.3733 1631.0000 2058.3733 -550383.81 -275191.88 -0.28739016 

27 1902.0000 1902.0000 1902.0000 299696.00 122398.00 0 . 13090003E-01 

28 1190.0000 1190.0000 1190.0080 73770.188 36885.099 0 .98533008E-02 

29 878.00000 878.00000 878.00000 157733.63 78866.813 0 . 13973827E-0 1 

30 681.00000 681.00000 681.00000 -109922.56 -52211.281 -0 . 1150028aE-01 

31 681.00000 681.00000 681.00080 192731.19 71365.563 0 . 1571 9287 E-01 

32 1902.0000 1902.0000 1902.8000 0.80800000 0.08000000 0.00000000 

33 1190.0000 1190.0000 1190.0000 0.80000880 0.00000000 0.00000080 

39 878.00000 878.00000 878.00800 0.80000800 0.00000000 0.00000080 

35 833.99536 0 .22913910E-03 0.21890225E-03 0.15999997 0 .79999983E-01 0 . 99227250E-01 

36 837.99536 0 . 18659518E-03 0 . 17865693E-03-0 .90800021E-01-0 . 20000011E-01 0 .26209069E-01 

37 893.99536 0.15919868E-03 0 . 15355896E-03-0 . 19999981E-01-0 . 99999905E-02 0.26006632E-01 

38 898.99536 0 . 19578699E-03 0 . 18996875E-03-0 . 19999981E-01-0 . 99999905E-02 0 .22053093E-01 

39 852.99536 0 .25921008E-03 0 . 29781892E-03-0 .39999989E-01-0 . 19999992E-01 0 . 13391986E-01 

90 859,99536 0 .29926179E-03 0 . 29356999E-03-0 . 39999999E-01-0 . 20000 0 00 E-01 0 . 10821126E-01 

91 859.99536 0 .27676396E-03 0 .27093795E-03 0,00000000 0.00000000 0 .22857092E-01 

92 1616.0000 1616.0000 1616.0000 976983.13 238991.56 0 .22137206E-01 

93 1631.0000 1631.0000 1631.0000 -139883.63 -67991.813 -0 . 62029966 E-02 

99 1902.0000 1902.0000 1902.0000 -331209.99 -165602.99 -0 . 177 17805E-0 1 

95 1190.0000 1190.0000 1190.0000 -218319.31 -109157.13 -0 . 19362775E-01 

96 878.00000 878.00000 878.00000 -127025.99 -63512.969 -0 . 10850735E-0 1 

97 681.00000 681.00000 681.00000 -101612.81 -50806.906 -0 . Ill 90S35E-01 

98 681.00000 681.00000 681.00000 0.80000000 0.00000000 0.00800000 

99 1902.0000 1902.0000 1902.0000 0.00000000 0.00000000 0.00000000 

50 1190.0000 1190.0000 1190.0000 0.00000000 0.00000000 0.00000000 

51 878.00000 878.00000 878.00000 0.00000000 0.00000000 0.00000000 

52 817.66260 0 . 37339996E-09 0 . 38963962E-09 0.00088000 0.00000000 -0 . 30239236 E-01 

53 817.66260 0 . 16920538E-03 0 . 17932208E-03 0 .39999999E-01 0 . 200 000 0 OE-01-0 . 12509651E-0 1 

59 819.66260 0 . 19990081 E-03 0 . 19913132E-03 0 .20000000E-01 0 . 99999979E-02-0 . 2237181 OE-0 1 

55 822.66260 0 .17758766E-03 0 . 18907662E-03 0 . 39999966E-01 0 . 1 9999981 E-01-0 . 1 9696929E-0 1 

56 827.66260 0 .23058079E-03 0 . 29095871E-03 0 .79999983E-01 0 . 39999992E-01-0 . 16818088E-0 1 

57 836.66260 0 . 226 0 9299E-03 0 . 23839192E-03 0 . 60000002E-0 1 0 . 30000001 E-01-0 . 39275722E-01 


58 

892.66260 

0.56570399E-03 

0.60062999E- 

03-0.23999995 

-0.11999995 

-0.93557775E-01 

59 

1616.0000 

1616.0000 

1616.0000 

0.00080000 

0.00000000 

0.00000000 

60 

1631.0000 

1631.0000 

1631.0000 

220892.99 

110996.19 

0.10157526E-01 

61 

1902.0000 

1902.0000 

1902.0000 

-56801.113 

-28900.555 

-0.30385756E-02 

62 

1190.0000 

1190.0000 

1190.0000 

-58973.597 

-29986.773 

-0.38798389E-02 

63 

878.00000 

878.00000 

878.00000 

-39558.890 

-19779.918 

-0.33791710E-02 

69 

681.00000 

681.00000 

681.00000 

-87399.125 

-93699.563 

-0.96259535E-02 

65 

681.00080 

681.00000 

681.00000 

-106129.31 

-53062.156 

-0.11687700E-01 

66 

1902.0000 

1902.0000 . 

1902.0000 

0.00800000 

0.00000000 

0.00000000 

67 

1190.0000 

1190.0000 

1190.0000 

0.00000000 

0.00000000 

0.00000000 

68 

878 .00000 

878 .00000 

878.00000 

0 .00000000 

0 .00000000 

0.00000000 

69 

9.7123709 

9.7123709 

9.7123709 

209.93999 

109.71997 

0.33333530E-02 

70 

209.93999 

209.93999 

209.93999 

-9.9528265 

-2.2269128 

-0.15995961E-05 





TABLE XVII. - INPUT FOR SUPPLY TEST CASE (MAINST) 


100 COnnOH/TEMIHT/TPRt1(2«) 

200 COMMO)I/STEPIT/FREQQ,FEEQ 

J 0 0 COnnOM/SUPLY /IJSS ( «> l , HOUEHG , PSOURC , TSOURC , U’SSCMO <t 1 , ftLERK 

0 0 COntlOM/PCBARR/PClMAX . PC2MAX , PC3MAX , PCAMAX , PCIMIN , PC2MIN . PC3MIN , 

500 IPCAMIM 

600 COtIMOM/ITCOMV/ VS (70),E(70).DELTAV(70),VSAVE(70),XD(<*1, 

700 1VDOTSV(70).XDCMV(A) 

800 COnnOM/PIM/ AD , AR , G , HH , HR , HC , AlOH , AWR , AWC , A J , CV , GAMMA , CVM , WSM ,H , 

900 IRST < 6 ) , VH . VR , VCOLD , UTOTT ( 9 1 , UTOT , UTOTSV , 

1000 2 ED,TOLSS,TOL2,TOLl,TOLDY,ISS,DEET,H,NMAX,VOE,VOC 

1100 3,V10T1,TUC,TCH,R0D1, 

1200 COMMOM/POUER/XDD ( 9 1 , PIHST , AMPLIT , TIME .OMEGA , ALPHA , XDDCNV ( 9 1 ,PS1 

1300 COMMOH/SAVIT/ VCOHV (70), VDOT (70), VGUESS( 7 0 ) , NODUG 

1900 1,VWS(70),VHIM(70),VDOTT(70) 

1500 COMMOH/AMGLS/SINl ( 9 ) , C051 ( 9 ) , SlHl 1 ( 9 ) ,TORQ ( 9 ) , TORQT 

1600 COMMOH/PRTIT/ VEl ( 9 ) , VCl ( 9 ) i 

1700 COMMOH/TEMPSO/TEO.THO, TRIO, TR20,TR30.TC0LD0, TOO, TUR10,TWR20,TMR30 

18 00 COMMOH/CARLOD/GTR AH , GR , AIE , AIUHEL , AIVEH , RIJHEEL , AINERT 

1900 COMMOH/COHST/PIE , DEGR , RANK 

2000 COMMOM/SUITCII/ ISUPST , 1C ALC , MATRIX , IPRTOP . IBR YTH , IHPCN V , ISPSTP 

2100 COMMOM/CYDATA/NPTPCY , ITRMAX , NUMBCY , STROKE , CYCLPR , SPDRPM, POSDEO 

2200 COMMOH/OUTPR/PCMAX , PCMIH , PEMAX , PEMIN , PAVEMP , KIJORK 

2300 COMMOH/MOVIT/DELPSI , FRAC , VDELTA , REF , MPAS , POSDEG . TOLPCG 

2900 C 

2500 C ENGIHE GEOMETRY 

2600 C 

2700 AD=23.7615 

2800 AR=1.1290 

2900 RODL=10.0 

3000 VR=115.332 

3100 VOE=11.959 

3200 VOC=32.5939 

3300 VH=52.9058 

3900 VCOLD=29.2378 

3500 STROKE=9.0 

3600 ALPHA=90.0 

3700 C 

3800 C HEATER 

3900 C 

9000 HH=71.9061 

9100 A11H=926.9669 

9200 TWH=705.0 

9300 C 

9900 C COOLER 

9500 C 

9600 HC=51.8739 

9700 A1JC=809.1279 

9800 TUC=86.0 

9900 C 

5000 C REGEHERATOR 

5100 C 

5200 HR=168.9762 

5300 AUR=66870.839 

5900 CVM=192.1395 

5500 USM=.6523 

5600 C 

5700 C TEMPERATURE DISTRIBUTION 

5800 C 

5900 TEO=629.99 

6000 THO=632.7778 

6100 TR10=505.5556 

6200 TR20=360.0 

6300 TH30=219.99 

6900 TCOLDO=105.0 

6500 TCO=105.0 

6600 C 

6700 C COHSTAHTS 

6800 C 

6900 RAHK=273.0 

7000 DEGR=57.296 

7100 PIE=3.1916 

7200 R=9125.6 

7300 G=10017.0 

7900 AO=1.0 



TABLE XVII. - Continued 


7500 

C 


7600 

c 

FLOW BESrSTANCES BETWEEN VOLUMES 

7700 

c 


7800 


EST(1)=. 38010 

7900 


BST(2)=. 76020 

8000 


EST(3)=. 76020 

8100 


RST(9)=. 76020 

8200 


RST(5)=. 76020 

8300 


RST(6)=. 38010 

8<i00 

c 


8500 

c 

HYDROGEN CONSTANTS 

8600 

c 


8700 


CP=1292.13 , 

8800 


GAMMA=1.39 

8900 

c 


9000 

c 

MATRIX INPUT DATA 

9100 

c 


9200 


VDELTA=.01 

9300 


FRAC=1 .0 

O'tOO 


TOL1=.001 

9500 


TOL2=.01 

9600 


TOLSS=.0001 

9700 


N=70 

9800 


NTMAX=70 

9900 


MPAS=50 

10000 


TOLPCG=.l 

10100 

c 


10200 

c 

SWITCHS 

1 0300 

c 


10900 


ICALC=1 

10500 


NONENG=l 

10600 


IPRTOP=200 

10700 


IBRYTH=1 

1 0800 


II!FCNV=0 

10900 


NOBUG=l 

11000 

c 


11100 

c 

RUN CONDITIONS 

11200 

c 


11300 


PSOURC=10.0 

11900 


TSOURC=10.0 

11500 


ALEAK=0.0 

11600 


SPDRPM=2900.0 

11700 


CyCLPR=5.00 

.11800 


SPDMAX=9000. 

11900 


CYCLPM=15.0 

12000 


POSDEO=270.0 

12100 


NCYSUP=10 

12200 


NCYSTP=200 

12300 

c 


12900 

c 

LOAD CONDITIONS 

12500 

c 


12600 


GTRAN=1.0/2.53 

12700 


GR=1.5 

12800 


WTENG-0.0 

12900 


WTI.'HEL=29.98 

13000 


WTCAR=1920.0 

13100 


RWIIEEL^30.98 

13200 

c 


13300 

c 

CYCLE DATA 

13900 

c 


13500 


NPTPCY-200 

1 36 0 0 


NUMBCY=10C 

13700 

c 


13800 

c 

CALCULATED INPUT 

13900 

c 


19000 


IF (KONENG .EQ. 0) N=68 

19100 


IF (IBRYTH .EQ. 1) TOLPCG=0.0 

19200 


REF=(TOL1+TOL2)/2.0 

19300 


ITRMAX=HPTPCY»'NUMDCY+1 

19900 


ISUrST=NCYSUP*'NPTrCY+l 

19500 


ISrSTP-NCYSTP'‘NPTPCY+l 

19600 


RANI'.=nA)!i;*96 0,0/273.0 



TABLE XVII. - Concluded 


l<t700 
14S00 
1<(900 
15000 
15100 
15200 
15300 
15900 
15500 
15600 
15700 
15800 
15900 
16000 
16100 
16200 
16300 
16900 
16500 
16600 
16700 
16800 
16900 
17000 
17100 
17200 
17300 
17900 
17500 
17600 
17700 
17800 
17900 
18000 
18100 
18200 
18300 
189 0 0 
18500 
18600 
18700 
18800 
1'8900 
1 9000 
19100 
19200 
19300 
19900 
19500 
19600 
1 9700 
19800 
19900 
20000 
20100 
20200 
20300 
20900 
20500 
20600 


AIE=0.0 

UTimEL=M™iEL/9.0 
RWIIEI,n=IU!HEEL/ 1 0 0 . 0 

AIWI1EL= ( UTU’HEL/2 . 0 1 « ( RMHELM»**2+ ( RUHELM/2 . 0 1 

AIVEH=WTC8R'«RUHELM'»*2 

fiIWHEL=WlJlIEi*60.56/.9279 

AIVEII=RIVEH*»1166 .56/131.9223 

RU11EEL=RMIIEEL«1 . 0/30 . 98 

CmH=2.59 

CMINS=Ct11H'<2.59 

CHIKC=CMI)IS»>2.59 

AD=RD/CMINS 

AR=AR/CniNS 

VR=VR/CniNC 

V0E=V0E/CM1HC 

V1I=VH/CHINC 

VCOLD=VCOLD/CMINC 

voc=voc/cmnc 

RODL=RODL/Ct1IN 

STROKE=STROKE/CniN 

TWl!=1.8^(TllH + 90.)-90. 

TllC=1.8’‘(TUC + 90.)-90. 

RIJlI=MJII/Ct1INS 
AUC=AWC/CniHS 
RUR = AI-!R/CHINS 
WSn=lJSn/.9536 

786/71. 9061 
HR = HR’»1. 861/168. 9762 
HC=HC«. 571/51. 8739 
TEO=1.8*(TEO+90.)-90. 

THO=1.8*(THO+90.)-90. 

TR10=1.8»*(TR10+9 0. 1-90. 

TR2O=1.8'‘(TR2O+90. 1-90. 

TH30=1.8*(TR30+90. 1-90. 

TCOLDO=l . 8^(TCOLDO+90 .1-90. 

TCO=1.8*(TCO+90. )-90. 

AJ=AJ«9337.92 
G=G’*386. 9/10017 . 

R=R'‘9197./9125.6 
RST(1)=RST( 1)^25./. 38010 
RST(6)=RST(6)»'25./.38010 
DO 1 1=2,5 

RST(I)=RST(I)*50./. 76020 
1 COMTIHUE 

CP=CP’*3. 9838/1292. 13 
CV=CP/GAnnA 

PSOURC=rSOURC'<2175 . /15 . 

TSOURC=l . 8«(TSOURC+90 . 1-90 . 

AiEAK=ALEAK« . 5/3 . 2258 
CYCLPR = CYCLrR'<2 1 7 5 . / 1 5 . 

CYCLPM=C YCLPM* 2 1 7 5 . / 1 5 . 

cvn=cvn**. 11/192. 1395 

T11R10=TR10 

T11R20=TR20 

TU'R30=TR30 

ISS=0 

MATRIX=1 

CALL ICSTIR 

STOP 

E1!D 



TABLE XVIII. - TEST CASE PRINTOUT 


STIRLING ENGINE FOUR CYLINDER. SEVEN VOLUMES PER CYLINDER, CONTROLS MODEL 


RUN CONDITIONS FOR THIS TRANSIENT 


TUII = 97J 

1.3 

TUC = 

359.3 

CYCLPR = 

5.000 

NDNENG 

= 

1 NUMBCY = 


100 


ALEAK = 0. 

0000 

PSOURC 

= 10.00 

TSOURC = 

283.3 

1SUP5T 

* 

2001 ISP5TP = 

00001 


TIME 

POMERT 

PFRICT 

PAUX 

PRRF 

PDRAC 

PIOAD 

PENG 

PNET 

ITRAN 


PAVEMP 


POSDEG 

TORQT 

TFRICT 

TAUX 

TRRF 

TDRAG 

TLOAD 

TENG 

TNET 

KUORK 


SPDRPM 

VKPH 

0.0000 

0.0000 

3.809 

2.627 

3.036 

0.078 

7.113 

6.075 

-13.59 


1 

0.028 


270.0 

50.00 

15.32 

10.05 

12.09 

16.22 

28.31 

25.77 

-0.8228E-01 


1 

2000. 

72.65 

0.2506E-01 

13.38 

3.808 

2.627 

3.036 

0.077 

7.113 

6.070 

-0.2111 


201 

5.060 


270.0 

52.68 

15.31 

10.05 

12.09 

16.22 

28.31 

25.76 

-1.391 


201 

2000. 

72.65 

0.5005E-01 

13.97 

3.803 

2.627 

3.036 

0.078 

7.113 

6.070 

0.3807 


001 

5.008 


270.0 

53.67 

15.29 

10.05 

12.09 

16.22 

28.31 

25.75 

-0.3830 


201 

2000. 

72.65 

0.7506E-01 

10.13 

3.803 

2.627 

3.036 

0.078 

7.110 

6.069 

0.5085 


601 

0.992 


270.0 , 

53.97 

15.29 

10.05 

12.09 

16.22 

28.31 

25.70 

-0.8202E-01 


201 

2000. 

72.65 

0.1001 

10.18 

3.801 

2.627 

3.036 

0.078 

7.110 

6.068 

0.6011 


801 

0.987 


270.0 

50.06 

15.29 

10.05 

12.09 

16.23 

28.31 

25.70 

0.1102E-01 


201 

2000. 

72.66 

0.1251 

10.21 

3.807 

2.627 

3.036 

0.078 

7.110 

6.073 

0.6196 


1001 

0.987 


270.0 

50.10 

15.31 

10.05 

12.09 

16.23 

28.31 

25.76 

0.6772E-01 


201 

2000. 

72.66 

0.1501 

10.22 

3.800 

2.627 

3.036 

0.079 

7.115 

6.071 

0.6307 


1201 

0.987 


270.0 

50.10 

15.30 

10.05 

12.09 

16.23 

28.31 

25.75 

0.7910E-01 


201 

2000. 

72.66 

0.1750 

10.22 

3.803 

2.627 

3.036 

0.079 

7.115 

6.070 

0.6392 


1001 

0.988 


270.0 

50.21 

15.29 

10.05 

12.09 

16.23 

28.31 

25.70 

0.1560 


201 

2000. 

72.66 

0.2000 

10.23 

3.803 

2.627 

3.036 

0.079 

7.115 

6.070 

0.6007 


1601 

0.989 


270.0 

50.21 

15.29 

10.05 

12.09 

16.23 

28.31 

25.70 

0.1513 


201 

2000. 

72.66 

0.2250 

10.20 

3.806 

2.627 

3.036 

0.079 

7.116 

6.073 

0.6080 


1801 

0.990 


270.0 

50.22 

15.30 

10.05 

12.09 

16.23 

28.32 

25.76 

0.1058 


201 

2000. 

72.66 



TABLE XVIII 


Continued 


0.2500 

19.29 

3.895 

2.627 

3.036 

9.080 

7.116 

5.973 

0.5519 

2001 

9.991 


270.0 

59.28 

15.30 

10.95 

12.09 

16.23 

28.32 

25.75 

0.2066 

201 

2900 . 

72.55 

0.2750 

19.78 

3.858 

2.627 

3.036 

9.080 

7.117 

6.986 

1.177 

2201 

5.028 


270.0 

55.28 

15.35 

10.95 

12.09 

16.23 

28.32 

25.80 

1.156 

201 

2900 . 

72.57 

0.3000 

15.09 

3.889 

2.628 

3.037 

9.081 

7.117 

6.517 

1.958 

2901 

5.112 


270.0 

55.17 

15.97 

10.95 

12.09 

16.23 

28.32 

25.93 

1.922 

201 

2901 . 

72.57 

0.3250 

15.39 

3.922 

2.628 

3.037 

9.082 

7.118 

6.550 

1.671 

2501 

5.199 


270.0 

55.99 

15.60 

10.95 

12.09 

16.23 

28.32 

26.06 

2.619 

201 

2901 . 

72.68 

0 . 3<*99 

15.55 

3.958 

2.628 

3.037 

9.082 

7.119 

6.586 

1.858 

2801 

5.285 


270.0 

57.82 

15.75 

10.95 

12.09 

16.29 

28.32 

26.20 

3.295 

201 

2901 . 

72.68 

0.3799 

15.78 

3.992 

2.629 

3.037 

9.083 

7.121 

5.621 

2.035 

3001 

5.370 


270.0 

58.52 

15.88 

10.95 

12.09 

16.29 

28.33 

26.33 

3.959 

201 

2901 . 

72.69 

0.3999 

15.98 

9.023 

2.629 

3.037 

9.089 

7.122 

6.652 

2.206 

3201 

5.953 


270.0 

59.91 

16.00 

10.96 

12.09 

15.29 

28.33 

26 . 96 ' 

9.628 

201 

2901 . 

72.69 

0.9299 

15.18 

9.058 

2.629 

3.038 

9.086 

7.123 

6.688 

2.366 

3901 

5.539 


270.0 

60.20 

16.19 

10.96 

12.09 

16.25 

28.33 

26.60 

5.268 

201 

2902 . 

72.70 

0.9999 

15.37 

9.086 

2.630 

3.038 

9.087 

7.125 

5.716 

2.527 

3601 

5.612 


270.0 

60.99 

16.25 

10.96 

12.09 

16.25 

28.33 

26.71 

5.895 

201 

2902 . 

72.71 

0.9798 

15.55 

9.116 

2.630 

3.038 

9.088 

7.127 

6.796 

2.682 

3801 

5.688 


270.0 

61.65 

16.36 

10.96 

12.09 

16.25 

28.39 

26.82 

6.990 

201 

2902 . 

72.71 

0.9998 

15.79 

9.199 

2.630 

3.039 

9.089 

7.128 

6.779 

2.833 

9001 

5.763 


270.0 

62.33 

16.97 

10.96 

12.09 

16 .26 

28.39 

25.93 

7.061 

201 

2902 . 

72.72 

0.5298 

15.91 

9.171 

2.631 

3.039 

9.091 

7.130 

6.802 

2.979 

9201 

5.839 


270.0 

63.00 

16.58 

10.96 

12.09 

16.26 

28.35 

27.09 

7.613 

201 

2903 . 

72.73 

0.5997 

17.08 

9.200 

2.631 

3.039 

9.092 

7.132 

6.831 

3.119 

9901 

5.909 


270.0 

63.68 

16.69 

10.96 

12.09 

16.26 

28.35 

27.15 

8.175 

201 

2903 . 

72.79 

0.5797 

17.29 

9.228 

2.632 

3.090 

9.099 

7.139 

6.860 

3.296 

9601 

5.972 


270.0 

69.32 

16.80 

10.96 

12.09 

15.27 

28.35 

27.25 

8.701 

201 

2903 . 

72.75 

0.5997 

17.90 

9.253 

2.632 

3.090 

9.095 

7.136 

6.885 

3.378 

9801 

6.038 


270.0 

69.95 

16.90 

10.96 

12.09 

16.27 

28.36 

27.36 

9.232 

201 

2909 . 

72.75 



TABLE XVIII. - Continued 


0.6246 

17.56 

4.278 

2.633 

3.041 

4.097 

7.138 

6.911 

3.508 

5001 

6.103 


270.0 

65.52 

17.00 

10.46 

12.09 

16.28 

28.36 

27.46 

9.699 

201 

2404 . 

72.77 

0.6496 

17.71 

4.303 

2.633 

3.041 

4.099 

7.140 

6.936 

3.630 

5201 

6.165 


270.0 

66.10 

17.09 

10.46 

12.09 

16.28 

28.37 

27.55 

10.18 

201 

2404 . 

72.78 

0.6745 

17.85 

4.329 

2.634 

3.042 

4.101 

7.143 

6.963 

3.748 

5401 

6.226 


270.0 

66.71 

17.19 

10.46 

12.09 

16.29 

28.37 

27.66 

10.69 

201 

2405 . 

72.79 

0.6995 

17.99 

4.353 

2.634 

3.042 

4.103 

7.145 

6.987 

3.861 

5601 

6.285 


270.0 

67.30 

17.28 

10.46 

12.09 

16.29 

28.38 

27.75 

11.17 

201 

2405 . 

72.80 

0.7244 

18.13 

4.378 

2.635 

3.043 

4.105 

7.147 

7.013 

3.974 

5801 

6.342 


270.0 

67.82 

17.38 

10.46 

12.09 

16.30 1 

28.38 

27.85 

11.59 

201 

2405 . 

72.81 

0.7493 

18.26 

4.396 

2.636 

3.043 

4.107 

7.150 

7.032 

4.081 

6001 

6.398 


270.0 

68.33 

17.45 

10.46 

12.09 

16.30 

28.39 

27.92 

12.03 

201 

2406 . 

72.83 

0.7743 

18.39 

4.418 

^.636 

3.044 

4.109 

7.152 

7.054 

4.188 

6201 

6.452 


270 .0 

68.83 

17.54 

10.46 

12.09 

16.31 

28.39 

28.00 

12.43 

201 

2406 . 

72.84 

0.7992 

18.52 

4.445 

2.637 

3.044 

4.111 

7.155 

7.082 

4.284 

6401 

6.505 


270.0 

69.32 

17.64 

10.46 

12.09 

16.31 

28.40 

28.10 

12.81 

201 

2407 . 

72.85 

0.8241 

18.64 

4.461 

2.638 

3.045 

4.113 

7.158 

7.098 

4.388 

6601 

6.556 


270.0 

69.82 

17.70 

10.47 

12.09 

16.32 

28.40 

28.17 

13.25 

201 

2407 . 

72.86 

0 .8491 

18.76 

4.481 

2.638 

3.045 

4.115 

7.161 

7.119 

4.482 

6801 

6.605 


270.0 

70.28 

17.78 

10.47 

12.09 

16.32 

28.41 

28.24 

13.62 

201 

2407 . 

72.88 

0.8740 

18.88 

4.504 

2.639 

3.046 

4.118 

7.163 

7.143 

4.570 

7001 

6.654 


270.0 

70.69 

17.86 

10.47 

12.09 

16.33 

28.42 

28.33 

13.94 

201 

2408 . 

72.89 

0.8989 

18.99 

4.519 

2.640 

3.046 

4.120 

7.166 

7.159 

4.667 

7201 

6.701 


270.0 

71.10 

17.92 

10.47 

12.09 

16.34 

28.42 

28.39 

14.29 

201 

2408 . 

72.90 

0.9238 

19.10 

4.538 

2.640 

3.047 

4.122 

7.169 

7.179 

4.750 

7401 

6.747 


270.0 

71.54 

17.99 

10.47 

12.09 

16.34 

28.43 

28.46 

14.64 

201 

2409 . 

72.92 

0.9487 

19.20 

4.556 

2.641 

3.047 

4.125 

7.172 

7.198 

4.832 

7601 

6.791 


270 .0 

71.90 

18.06 

10.47 

12.09 

16.35 

28.44 

28.53 

14.95 

201 

2409 . 

72.93 

0.9736 

19.30 

4.573 

2.642 

3.048 

4.127 

7.175 

7.215 

4.914 

7801 

6.835 


270.0 

72.32 

18.13 

10.47 

12.09 

16.36 

28.44 

28.60 

15.28 

201 

2410 . 

72.95 



TABLE XVIII. - Continued 


0 . 99 R 5 

19.90 

9.592 

2.693 

3.099 

9.130 

7.178 

7.235 

9.990 

8001 

6.876 


270.0 

72.72 

18.20 

10.97 

12.09 

16.36 

28.95 

28.67 

15.60 

201 

2910 . 

72.96 

1.023 

19.50 

9.608 

2.693 

3.099 

9.132 

7.182 

7.252 

5.063 

8201 

6.917 


270.0 

73.08 

18.26 

10.97 

12.09 

16.37 

28.96 

28.73 

15.90 

201 

2911 . 

72.98 

1.098 

19.59 

9.625 

2.699 

3.050 

9.135 

7.185 

7.269 

5.191 

8901 

6.957 


270.0 

73.93 

18.32 

10.97 

12.09 

16.38 

28.96 

28.79 

16.18 

201 

2911 . 

72.99 

1.073 

19.68 

9.692 

2.695 

3.051 

9.137 

7.188 

7.287 

5.205 

8601 

6.995 


270.0 

73.78 

18.38 

10.97 

12.09 

16.38 

28.97 

28.86 

16.96 

201 

2912 . 

73.01 

1.098 

19.77 

9.659 

2.696 

3.051 

9.190 

7.101 

7.305 

5.273 

8801 

7.033 


270.0 

79.13 

18.95 

10.98 

12.09 

16.39 

28.98 

28.92 

16.73 

201 

2912 . 

73.02 

1.122 

19.85 

9.671 

2.697 

3.052 

9.193 

7.195 

7.317 

5.339 

9001 

7.069 


270.0 

79 .99 

18.99 

10.98 

12.09 

16.90 

28.98 

28.96 

17.00 

201 

2913 . 

73.09 

1.197 

19.93 

9.686 

2.697 

3.053 

9.195 

7.198 

7.333 

5.901 

9201 

7.109 


270.0 

79.79 

18.59 

10.98 

12.09 

16.90 

28.99 

29.02 

17.28 

201 

2913 . 

73.05 

1.172 

20.02 

9.700 

2.698 

3.053 

9.198 

7.201 

7.399 

5.965 

9901 

7.138 


270.0 

75.11 

18.60 

10.98 

12.09 

16.91 

28.50 

29.08 

17.59 

201 

2919 . 

73.07 

1.197 

20.09 

9.715 

2.699 

3.059 

9.151 

7.205 

7.369 

5.523 

9601 

7.171 


270.0 

75.90 

18.65 

10.98 

12.09 

16.92 

28.50 

29.13 

17.76 

201 

2919 . 

73.09 

1.222 

20.17 

9.729 

2.650 

3.055 

9.159 

7.208 

7.379 

5.587 

9801 

7.209 


270.0 

75.60 

18.68 

10.98 

12.09 

16.93 

28.51 

29.16 

17.92 

201 

2915 . 

73.10 

1.296 

20.29 

9.791 

2.651 

3.055 

9.157 

7.212 

7.592 

5.639 

10001 

7.236 


270.0 

75.88 

18.75 

10.98 

12.09 

16 .93 

28.52 

29.23 

18.13 

201 

2915 . 

73.12 

1.271 

20.31 

9.759 

2.652 
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Figure 1. - Four working space model. Last number of a variable name indicates working space; cylinders are numbered in the order that they reach top dead center; primed variables Indicate intervolume temperatures, 
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Figure 2. - Stirling engine drive dynamics 




































Figure 4. - Overall simulation structure. 
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Figure 5. - Simulation results; comparison of complex (seven vol- 
ume) and simplified (three volume) simulations. 
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Figure 8. - Concluded. 





1. Report No. 2. Government Accession No. 

NASA TM -83053 

3. Recipient's Catalog No. 

4 . Title end Subtitle 

A Four -Cylinder Stirling Engine Computer Program 
With Dynamic Energy Equations 

5. Report Date 

May 1983 

6. Performing Organization Code 

778-35-03 

7. Authorlsl 

Carl J. Daniele and Carl F. Lorenzo 

8. Performing Orgenitation Report No. 

E-1348 

10. Work Unit No. 

9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

12. Sponsoring Agency Name and Address 

U. S. Department of Energy 
Office of Vehicle and Engine R&D 
Washington, D. C. 20546 

14. Sponsoring Agency ftx^Report NO. 

DOE/NASA/51040-44 

IS. Supplementarv Notes 


Final report. Prepared under Interagency Agreement DE-AIOl -77CS51040. 


16. Abstract 

A computer program for simulating the steady -state and transient performance of a four- 
cylinder Stirling engine is presented. The thermodynamic model includes both continuity and 
energy equations and linear momentum terms (flow resistance) . Each working space betweeir 
the pistons is broken up into seven control volumes. Drive dynamics and vehicle load effects 
are included. The model contains 70 state variables. Also included in the model are piston- 
rod-seal leakage effects. The computer program includes a model of a hydrogen supply 
system, from which hydrogen may be added to the system to accelerate the engine. Flow 
charts are provided. 


17. Key Words (Suggested by Author(sM 

Stirling engine 

Controls 

Modeling 

Computer program 

18. Oistribution Statement 

Unclassified - unlimited 
STAR Category 34 
DOE Category UC-96 

19. Security Oassif. (of this report) 20. Security Clasiif. (of this page) 

Unclassified Unclassified 

21. No. of Pages 

22, Pilet* 

AOS 


‘ For sale by the National Technical Information Service. Springfield, Virginia 22161 
















End of Document 



